A Cellular Automata Model for Citrus Variagated Chlorosis 
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A cellular automata model is proposed to analyze the progress of Citrus Variegated Chlorosis epi- 
demics in Sao Paulo oranges plantation. In this model epidemiological and environmental features, 
such as motility of sharpshooter vectors which perform Levy flights, hydric and nutritional level of 
plant stress and seasonal climatic effects, are included. The observed epidemics data were quantita- 
tively reproduced by the proposed model varying the parameters controlling vectors motility, plant 
stress and initial population of diseased plants. 
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I. INTRODUCTION 

Citrus Variegated Chlorosis (CVC) is an economically 
relevant disease affecting citrus Q . In the Sao Paulo re- 
gion (Brazil), one of the important citrus growing areas 
of the world, responsible for about 30% of the world pro- 
duction, CVC reduces the size and number of fruits by 
more than 35% CVC is considered to be potentially 
the most devastating citrus disease and represents the 
main threat to the Brazilian citrus industry, with annual 
revenues of the order of 1.2 to 1.4 billions of dollars. The 
losses associated with the disease are estimated in about 
100 millions dollar yearly Q. 

CVC is caused by a xylem-limited bacterium, Xyl- 
lela fastidiosa transmitted by xylem feeding, suc- 
torial sharpshooter leafhoppers (Hemiptera: Cicadelli- 
dae) 1^,^. In Sao Paulo, the species Dilohopterus costal- 
imai appears to be the most efficient vector for CVC 
transmission Nowadays, a sweet orange cultivar re- 
sistant to X. fastidiosa is unknown, and control practices 
for CVC (bactericidal agents, systematic pruning of in- 
fected branches, chemical control of vectors, and/or roug- 
ing of severe affected plants) are expensive, ineffective, or 
environmentally damaging. 

Recent studies on various aspects of the epidemiology 
of CVC ( 01 , and references therein) have provided fun- 
damental information which can be used to develop a 
cellular automata (CA) model of the pathosystem. CA 
or other epidemic models could become relevant tools to 
address numerous practical and experimental questions: 
forecasting the progress and final intensity of CVC, plan- 
ning and evaluation of strategies for disease control and 
determination of the relevant mechanisms involved in the 
disease spreading. 

In this paper, we propose a simple CA model to simu- 
late the CVC progress in which some epidemiological and 
environmental features, such as vectors motility, plant 



stress and seasonal modulations are included. The sim- 
ulational results are compared with the CVC progress 
curves in time and spatial infection patterns observed in 
the Sao Paulo region. 



II. EXPERIMENTAL DATA ON CVC EPIDEMICS 
A. CVC progress in time 

The CVC epidemic was observed by visual assessments 
of typical symptoms occurring on leaves or fruits, in 
eleven groves of Pera, Hamlin and Natal sweet oranges 
cultivated in two farms of the northern areas (Bebedouro 
and Colina counties) of Sao Paulo state, Brazil. In 
such areas, severely attacked by CVC, are planted the 
more susceptible cultivars having the supposedly more 
propitious age for disease development. The field data 
were collected over a 20-month period, from Septem- 
ber, 1994 through March, 1996. The CVC incidence 
were bimonthly mapped and the data for each area and 
each evaluation were transformed to proportion of symp- 
tomatic plants for temporal characterization of the dis- 
ease spreading. 

The CVC progress curves are shown, for four differ- 
ent groves, in figure I. AU of them are double sigmoid, 
which is a clear indication that CVC is a polycyclic dis- 
ease characterized by the existence of two phases: one 
in which the disease spreading is fast, contrasted by an- 
other one in which the epidemic development is almost 
stopped. For each grove the observed data sets are fitted 
by five parameter logistic curves B of the form 



N{t) 



Pi 



1 + exp[-{p2 P3t + Pif^ +P5t^)] ' 



(1) 



Table 1 gives the corresponding parameters and the co- 
efficients of determination (i?^) have been listed. In ad- 
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dition to the coefficient, the residual sums of squares 
for error and the consistency of the predict values for the 
upper asymptotic fraction of diseased plants (pi < 1) 
were take into account to select the logistic among the 
Gompertz and monomolecular generalized models with 
four or five parameters. It is important to notice that 
sigmoid curves can be generated by distinct models (fit- 
ting equations). Indeed, the temporal increase of citrus 
tristeza virus, whose vectors are aphid species, follows a 
nonlinear Gompertz model 0. 



binary symbol (Ti{t), describing the plant state (o'i(t) = 0: 
normal; ai(t) = 1: diseased) is identificated with a step 
(to the right or to the left) of a one-dimensional walk. 

Specifically, to an unique spatial pattern 
{(Ti(t), (T2{t), . . . , (7Ar(t)} of N plants at fixed time t, cor- 
responds a profile given by the sequence of the walker dis- 
placements hi after i unit steps, {/ii(i), h2{t), . . . ,h]\[{t)}, 
defined as 



(2) 



B. Infected cluster size distribution functions. 

On any one of the orangeries containing about lO'^ trees 
there is no unique inoculum source. Each single infected 
plant or small initial group of infected plants grows by 
inoculating its adjacent neighbors, and aggregates with 
other affected trees, forming large clusters. As a result, 
the mean cluster size of infected trees increases in time 
and, in order to describe the disease spreading, it is nec- 
essary to investigate the dynamic aspects of the distri- 
bution of infected plant aggregates generated by CVC 
progress. 

A cluster of diseased plants have been defined as any 
set of interconnected infected trees which are spatially 
isolated from any other group of diseased plants in the or- 
angery. Then, the cluster size distribution function (t) 
is the fraction of these clusters consisting of s infected 
plants at time t, and were directly obtained by counting 
the number and size of diseased clusters present in the 
spatial patterns of CVC epidemic, like those shown in 
figure H, at each observation time. 

Figure ^ shows the distribution Ug {t) for one of the ob- 
served orangeries. It suggests that ns{t) follows a power- 
law distribution, that is, ns(i) ^ s over at least one 
decade of their argument s at any given time of CVC 
progress. A power-law distribution indicates the absence 
of a characteristic scale for the size of diseased plant ag- 
gregates. The exponents describing the power-law de- 
cay of ns(t) have values between 1.4 and 1.8 for all the 
studied or ang eries. These values are characteristic of an 
1 // noise , as are commonly referred the 1/ spectra 
with 7 in the interval [0, 2]. Therefore, the CVC infection 
dynamics has an l//-like signal, which, from the physi- 
cal point of view, is a sign of a cooperative phenomenon 
occurring in a spatially extended nonequilibrium system. 



where pj = 1 if ctj = 1 (step to the right), or pj = — 1 
if (jj — (step to the left). Profiles generated at two 
distinct observation times in a given orangery using this 
walk process are shown in figure ^. 

Obtained the profiles by the walk process, we can inves- 
tigate the nature of their correlations through the anal- 
ysis of the profile roughness |^ . The statistical measure 
W, which characterizes the roughness of the walk profile, 
is defined by the rms fluctuation in the displacement 



W{N,t) = 



1 ^ 



where 



N 



(3) 



(4) 



is the mean displacement of the walk. 

For self-affine profiles the roughness W{N) will be de- 
scribed by a power law scaling 



W{N) - N 



H 



(5) 



with the exponent H restricted to the interval [0, 1] 
and related to the fractal dimension of the profile 
H = 1/2 corresponds to a random walk; H > 1/2 im- 
plies that the profile presents persistent correlations, and 
profiles with H < 1/2 are anticorrelated. 

As one can see in figure ^, the profiles generated from 
the CVC epidemic patterns usually have drifts. This is 
the reason why we use the method of roughness around 
the rms straight line [l^ to evaluate the Hurst exponent 
H. In this method the roughness W{N, e) in the scale e 
is given by 



C. Self-afHne profiles in CVC evolution patterns. 



1 ^ 



(6) 



Self-affine profiles g,|9| can be generated from the spa- 
tial patterns of diseased plants using various methods. 
The simplest of them is the 1:1 mapping between a given 
spatial configuration at time t^such as those shown in fig- 
ure ^, and a "walk process" [|lO|,0 . In this method each 



and the local roughness ^^(e) is defined as 



\ 
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2e+l 



j=i-<i 



-h{e)]V. (7) 
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ai{e) and bi{e) are the linear fitting coefficients to the 
displacement data on the interval [i — e, j + e] centered 
on the site i. Again, self-afRne profiles satisfy the scaling 
law 

T^(e)~e^ (8) 

The method described here was used to characterize 
the spatio-temporal patterns generated by elementary 
one-dimensional deterministic cellular automata [ p^ . 

A Typical log-log plot of W{e) x e used to calculate 
the Hurst exponent is shown in figure ^ This exponent 
characterizes the spatial patterns generated by CVC epi- 
demic. For all the analyzed orangeries the profiles are 
self-affine with Hurst exponents markedly distinct from 
1/2, which means long-range correlations present in the 
spatial disease patterns. Also, the roughness exponents 
increases from a initial value around 1/2, indicating a 
random infection pattern at the beginning, towards its 
maximum value 1 which corresponds to a totally infected 
orangery. Thus our results show that CVC epidemic gives 
rise to aggregated patterns in which the inoculum level 
tends to be high in scattered sequences of neighboring 
plants, i. e., infected plants tends to be close to an- 
other diseased trees and the same holds for normal plants. 
Therefore, it appears that the pathogen occurs in small 
clusters which progressively expands by a contagion pro- 
cess mediated by vectors which predominantly spreads 
from plant to plant. 



III. A CELLULAR AUTOMATA MODEL FOR 
CVC 

Stochastic CA models were used before in plant pathol- 
ogy to simulate diseases spreading through spore disper- 
sal p6| , infection dynamics of R. Solan i |ll7|l a nd infection 
of cereal roots by the take-all fungus [|18| , |19|| , but tradi- 
tionally the mathematical modeling of plant disease is 
based on systems of differential equations |l^ . 

CAs are totally discrete dynamical systems (discrete 
space, discrete time and discrete number of states) which 
provide simple models for a great number of problems 
m science ||l|Jl|. With each site (noted i), IS associ- 
ated a variable ai, which can be in K different states 
ai = 0,1, K — 1. The dynamics is defined, at each 
time step, by rules depending on the values at previous 
time of {ci} associated with a given number of q arbi- 
trary sites (called inputs). Usually one considers regu- 
lar lattices and the inputs refer to the sites on the local 
neighborhood only. The local rules of a CA may be prob- 
abilistic or deterministic and the sites are simultaneously 
updated. 

In order to design the CA model for CVC spreading, 
we take into account the following basic features of CVC 
pathosystem characterized in the previous section. The 
bacterium X. fastidiosa is transmitted by sharpshooter 



vectors of rather limited motility in the groves. This 
hypothesis is coherent with the results obtained for the 
roughness or Hurst exponents describing the CVC infec- 
tion patterns. In diseased plants, the bacterium is sys- 
temic, nutritional unbalance and general weakness are 
commonly observed. Thus, the infected sites continu- 
ously act as inoculum sources to other healthy plants. 
Since there was no measured effect of wind direction or 
machine based cultural practices on CVC spreading [Q, 
the vectors flies were assumed to be completely random 
in our model. Also, in healthy stressed trees the ob- 
served sharpshooter population is small, since these in- 
sects are preferable attracted by plants with new vege- 
tative growth. Finally, as shown by the CVC progress 
curves, seasonal effects play a central role in disease 
spreading. In fact, the fastest CVC spreading progress 
is observed from September/October (flowerage) through 
March (end of summer) , a period associated to high tem- 
peratures, regular rains and vegetative growth. The sea- 
sonal modulations are included in the CA model through 
the variation in the motility of sharpshooter vectors as 
well as in the fraction of normal plants under hydric and 
nutritional stress. The functional form assumed to model 
such seasonal variations in our model is a sine wave func- 
tion. 

In our CA model the orangery is represented by square 
lattices of linear size L with null fixed boundary condi- 
tions (isolated grove, i. e., all the state variables are 
zero outside the lattice). The state of each plant is de- 
scribed by two binary variables: a'l^j and cr^j^ where 
i,j — 1,2,..., L. a^^J = represents a healthy (normal) 
plant and til^^ = 1 a diseased plant. On the other hand, 

(2) 

al ■ — 1 represents a plant under hydric or nutritional 

(2) 

stress and ai ■ = a non-stressed plant. Finally, the frac- 
tion of inoculative vectors in each site is also represented 
by a binary variable Ti^. Tij — means a low fraction of 
inoculative sharpshooter leafhoppers in the insect popu- 
lation and Tij = 1 a high inoculative fraction. 

In all simulations, random initial conditions in which 
any plant of the lattice is diseased, with probability , 
and stressed, with probability Pst(O), have been used. 
Since CVC causes severe stress in a diseased plant [0, 

a infected tree (a^^j — 1) becomes immediately stressed 

(dj y = 1) in our model. 

All the sites are simultaneously updated using the fol- 
lowing local rules: 

i) For an infected site the correspondent state at the 
next time step is a^^j = 1, af''- — 1 (diseased and stressed 
tree) and Ti_j — 1 (high fraction of inoculative vectors). 
In addition, as shown in figure 0, each infected site acts 
as an inoculum source for n„(fy distinct plants at dis- 
tance rk, k — 1,2, . . . ,nv, chosen at random accordingly 
a symmetric Levy distribution 
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p{r) = ^J dt exp{it{pL - r) - It^). (9) 

Thus, the lengths of each one of the vector flights 

are not constant but rather are chosen from a proba- 
bility distribution with a power-law tail. Each one of 
these selected neighbors will assume, at the next time 
step, the state ct - = 1, u''-^^ = 1 and = 1, if it is a 
normal and non-stressed plant. Otherwise, the selected 
neighbor will stay in the same previous state. Thus, in 
our CA model a healthy stressed tree is not infected by 
CVC, since their main vectors (Dilobopterus Costalimai 
and Acrogonia sp.) are preferentially observed in plants 
exhibiting young buds and leaves. In contrast, a healthy 
and non-stressed plant becomes diseased if it is reached 
by inoculative vectors coming from at least one infected 
site. 

ii) For a normal (stressed or non-stressed) plant not 
target by a given diseased plant, the correspondent state 
at the next time step is a^^j = (normal) and Ti.j = 
(low fraction of inoculative vectors). Yet, the value of 

(2) 

CTj y = or 1 is chosen at random with probability 
Pst = 1 -Pnst{t). 

iii) In order to simulate the seasonal effects on both, 
plant stress and vectors motility, the number of inocula- 
tive vector flights n^(t), and the probabihty associated 
to a non-stressed plant Pnst{t) are periodic functions of 
time given by: 

n,{t)^no + INT[nostn{— + ^)]; (10) 



Pnstit) ^ Prmn + - '^Pmin) X ' , (H) 

where INT means the integer part, Uq is the time av- 
erage number of inoculative vector flights, Pmin is the 
minimum probability to find a non-stressed plant, and T 
is the period of one year. </) is a common phase angle to 
describe possible time shifts between the simulated and 
field data. 

At last, we shall discuss some simplifications of our 
CA model for CVC progress. The use of a symmetric 
periodic function to model the seasonal effects should 
be thought as a rough approximation to the much more 
complex climatic variations observed in nature. Another 
simplification is that once infected a given plant imme- 
diately acts as an inoculum source, in contrast to the 
classical notion of a discontinuous infectious period. It is 
known that the spread of disease involves the interplay 
of two dynamical processes: the mechanisms of transmis- 
sion and the evolution of the pathogen within hosts. The 
basic questions on how the bacterium spread within the 
xylem system and what is the mechanism of pathogenesis 
in CVC are unanswered. Particularly, it seems that the 
CVC symptoms in plants depend on the rate and extent 



of colonization by the bacteria. Also, a recent study pQ ] 
shows that an efficient transmission of Xylella fastidiosa 
by vectors occurs only after its population overcome a 
threshold in plant hosts. In addition, the transmission 
rate increases as the bacterial population in a plant in- 
creases. These features are not included in our CA model 
due to the lack of information about the system Xylella 
fastidiosa in citrus. 

IV. RESULTS 

Now we shall report on the simulational results for our 
CA model. In all simulations a linear size L — 200, a 
seasonal period T = 12 months (one year) and a phase 
angle (j> = —60° (a shift of two months between the ob- 
served and the simulation initial times) were used. The 
remaining five CA parameters, namely, no, Pm/, Pmim 
/i, and a were varied in order to compare the simulated 
and observed CVC progress curves. The results for four 
groves are shown in figure 0. As one can see, the quali- 
tative behavior of the CA progress curves show the same 
functional behavior of the measured curves. Even a sur- 
prisingly quantitative agreement between the CA and the 
observed CVC progress curves was obtained. The CA 
parameters values are listed in Table 2. Therefore our 
simulation results suggest that the Levy distribution of 
vectors flights is universal, with fi — 4 and a = 0.68. 
The same holds true for no, the time average number 
of sharpshooter vectors flies starting from each infected 
plant, which assumed a value 5 for all four simulated 
groves. 

At this point, it is interesting to note that our CA 
model predicts even triple, quadruple or superior sigmoid 
progress curves depending on the time elapsed up to a 
total infection of the grove. This simulational predic- 
tion can be easily tested simply by observing the CVC 
progress in field for a period greater than 20 months 
as done in the present work. Moreover, such triple 
or quadruple sigmoid curves could be mathematically 
modeled by using generalized logistic, monomolecular or 
Gompertz functions only if the disease progress curves 
were subdivided into three or four parts that are ana- 
lyzed separately. However, this approach is inadequate 
to describe the entire disease dynamics and to determine 
several parameters of epidemiological importance pl| , p2| . 

In figure ^, it is shown a simulated temporal sequence 
of CVC incidence maps qualitatively similar to the ob- 
served spatial patterns of CVC (see figure ^) . From such 
incidence maps, one can determine the dynamic infected 
cluster size distribution function, ns{t), and the rough- 
ness exponent H characterizing the self-affinity of CVC 
infection profiles. Figures || and |l^ show, respectively, 
the cluster size distribution functions ns(t), and typical 
log-log plots of W{e) X e corresponding to the simulated 
infection maps for the grove SJ71 at various observation 
times. Both, the power-law decay of the infected cluster 
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size distribution functions and the roughness exponents 
characterizing the spatial disease patterns, indicate the 
presence of long-range correlations in the CVC develop- 
ment. Roughness exponents greater than 1/2 mean that 
in the neighborhood of a diseased plant the probability 
to find another infected tree increases. A significative 
aggregation of diseased plants suggest that the pathogen 
predominantly spreads from plant to plant. Since our CA 
model permits fast simulations of many large samples of 
artificial CVC pathosystems under various epidemiolog- 
ical contexts, the numerical values for the exponents a 
and H can be determined with a reliable statistical pre- 
cision difficultly attained in actual field observations. 

It is important to emphasize that a simple random 
walk distribution for the inoculative vector flights con- 
strained to a local neighborhood of radius r(t) around 
each diseased plant, also generates progress curves in 
good agreement with the field data. However, the re- 
sulting CVC incidence maps are clearly distinct from the 
observed ones. Random flights of inoculative vectors pro- 
duce spots of diseased plants artiflcially isolated in space, 
which grow up to merge with other infected clusters. In 
contrast, a Levy distribution permits rare long range vec- 
tor flights which appears to be an essential feature to 
explain the scale invariance observed in CVC spreading. 
Indeed, conventional random walks used to model forag- 
ing behavior in biology ||2^ predict a Poisson instead of 
a scale- invariant power-law distribution. Thus, our re- 
sults suggest that the inoculative sharpshooter leafhop- 
pers (~ 1 cm in size) perform long flights of random for- 
aging, searching for non-stressed plants with new vegeta- 
tive growth unpredictably dispersed over several square 
kilometers. A possible explanation is that for insects 
operating in swarms or flocks comprised of N walkers, 
Levy flight search patterns (for which Nt sites are vis- 
ited after t steps) are much more efficient than Brownian 
walk foraging patterns (for which only tln{N /Int) dis- 
tinct sites are visited) [E5[. It is interesting to mention 
that Levy flight search strategies are also observed in al- 
batrosses [ p6[ . 

V. CONCLUSIONS 

In this study the spatio-temporal analysis of CVC 
spreading was carried out. The shape of the observed 
CVC progress curves was double sigmoid best fitted by 
five parameters generalized logistic function. This means 
that CVC is a polyciclic disease in which a phase of rapid 
progress alternates with another one of almost paralyzes. 
In addition, both the power-law decay of the infected 
cluster size distribution functions and the roughness ex- 
ponents characterizing the spatial disease patterns, indi- 
cate the presence of long-range correlations in the CVC 
development. 

In order to understand the basic mechanisms by which 
the features discussed above emerge, a CA model was 



proposed. It takes into account the motility of sharp- 
shooter vectors, the hydric and nutritional level of plant 
stress, as well as seasonal climatic effects. Varying the 
CA parameters controlling these factors, a good agree- 
ment among simulational and all the observational data 
was achieved, suggesting that the actual relevant mech- 
anisms of CVC spreading were really captured and evi- 
denced by the evolution rules of the proposed CA model. 
Therefore, our model suggest that the average number 
n.v(t) of Levy flights performed by the sharpshooter vec- 
tors as well as plant stress, described by the probability 
Pst, are the most fundamental parameters determining 
the aspects of CVC spreading. Also, all of them are af- 
fected by seasonal variations. 
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Table 1 Coefficients of determination and esti- 
mated values of the five parameters describing the gener- 
alized logistic functions used to fit the observed data of 
CVC progress in four of the studied groves. 



Parameters 


Pi 


P2 




Pi 


Ps 




SJOl 


0.223 


-6.000 


1.647 


-0.158 


0.005 


0.989 


S,J()7 


0.!)()8 


-G. ill 


i. ()(>:) 


-0.15:! 


0.005 


0.996 


SJ71 


1.000 


-4.854 


1.307 


-0.126 


0.004 


0.989 


SJ75 


0.526 


-6.638 


1.396 


-0.120 


0.004 


0.993 



Table 2 CA parameters used to simulate CVC 
progress curves corresponding to the observed data for 
four groves. 



Parameters 




a 


no 


Pin/ 


Pmin 


SJOl 


4 


0.68 
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0.0002 


0.001 


SJ67 


4 


0.68 


5 


0.04 





SJ71 
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0.68 
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0.06 


0.35 


SJ75 


4 


0.68 


5 


0.009 
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FIG. 1. Observed CVC progress curves in four groves SJOl, 
SJ67, SJ71 and SJ75. The curves represent generalized five 
parameter logistic fittings to the observed data. 
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FIG. 2. Temporal sequence of CVC incidence maps in the 
grove SJ 71 with 3960 Pera orange trees. Six representative 
evaluations (11/94, 01/95, 05/95, 09/95, 11/95 and 01/96) are 
shown. Each black square corresponds to one symptomatic 
plant. 
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FIG. 3. Observed cluster size distribution function ns{€) 
of clusters containing s diseased plants for the grove SJ71 at 
the observation time f = 6 months. The straight line repre- 
sent the best fit to the data; its slope gives the exponent a 
describing the power-law decay of ns{t). 
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FIG. 4. Self-affine profiles generated by the walk process 
method at two distinct observation times in the grove SJ75, 
containing 2,880 Pera orange plants. The rescaled profile at 
t = 14 months was obtained by subtracting the corresponding 
linear fitting of the data. The inset is an enlargement of the 
central region of the profile at f = 14. 
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FIG. 5. Typical log- log plot of W(e) x e used to calculate the roughness or Hurst exponent characterizing a rough profile, 
in this case the CVC spatial patterns in grove SJ75 at two distinct observation times. The Hurst exponent corresponds to the 
slope of the straight-line fitting the linear part of the curve. 




FIG. 6. Schematic representation of CVC spreading from a diseased plant (central gray square) used in our CA model. The 

central infected site act as an inoculum source for riv (t) = 4 other distinct sites whose distances r arc chosen from a symmetric 
Levy distribution. These "target" sites are reached through the vectors files represented by the arrows. n„ (t) change seasonally 
with time. 
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FIG. 7. Comparison between the CA generated and fitted CVC progress curves for groves (a) SJOl, (b) SJ67, (c) SJ71 and 
(d) SJ75. The observed data are also shown, and the used CA parameters are those listed in Table 2. The data correspond to 
an average over 100 diflferent realizations of CA evolution. The error bars are equal or smaller than the symbol syze. 
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FIG. 8. A simulated sequence of CVC incidence maps in the grove SJ71. Each black square corresponds to one symptomatic 
plant. 
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FIG. 9. Simulated cluster size distribution function ns{t) for the grove SJ71 at various observation times. The straight line 
represent the best fit to the data; its slope gives the exponent a describing 
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FIG. 10. Typical log- log plot of W[t)y.t corresponding to the simulated infection profiles for the grove SJ71 at two observation 
times. The straight line represent the best fit to the data; its slope gives the roughness exponent H describing the profile. The 
data correspond to an average over 20 different realizations of CA evolution. 
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